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£\j ' We have explored the magnetosphere described by a numerical solution for an axisym- 

metric rotating field satisfying the force-free and ideal MHD conditions everywhere. The 
electric current distribution is determined by the requirement of continuity and smoothness 
of the flux function across the light cylinder. The overall magnetic field structure is obtained 
using the numerical result. We also checked the drift velocity of a charged particle, which 
is caused by the resultant electric and magnetic fields. It is found that the velocity exceeds 
the light speed beyond a few times the size of the light cylinder. This fact suggests the 
breakdown of some underlying assumptions. 
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qq ■ §1. Introduction 

It is well known that rotating magnetized neutron stars emit radiation from radio 
waves to gamma rays, and constitute a power supply to surrounding nebula. The 
C<*) \ basic idea of a pulsar was proposed more than 30 years ago (e.g., Refs.tTJ and©), 

but a detailed consistent model has not yet been constructed. Almost all the energy 
r S^ ■ is carried out from a pulsar by electromagnetic fields, but it is converted to particle 

kinetic energy at large distances.^ Several authors have discussed the conversion 
mechanisms, such as non-ideal MHD effects,^ magnetic reconnectiorf^® and plasma 
-£3 . instability.^^ Because the dissipation of electro-magnetic fields is related to the 

acceleration of particles, and hence the observed radiation, the determination of the 
mechanism is very important. 

The global structure of the magnetosphere, even without dissipation, is not easily 
determined. The configuration is strongly coupled with the plasma flow. Several 
authors have used simplified electromagnetic models, such as a monopole magnetic 
field, to investigate the plasma dynamics. For instance, MicheP^'^3 considered cold 
MHD winds for an aligned pulsar, and Bogovalov^^ extended it to the oblique case. 
The monopole magnetic field itself is unrealistic, but some results may be used as a 
model of pulsar magnetosphere. It is desirable to consider the plasma dynamics with 
a more realistic magnetic field, e.g., a dipole field near the star. We may gradually 
confirm the existing general picture, or we may find incorrect aspects of previous 
models, using a realistic model. 

The pulsar equation determines an axially symmetric MHD field in the force- 
free limit (e.g., Refs. HI ]) .ITS ]) JT^ ll ) . Recently, Contopoulos et alP^ obtained a new 
type of numerical solution of the pulsar equation. We refer to their solution as 
the CKF solution in this paper. Their solution represents the interesting structure 
that the magnetic field is a dipole near the star and asymptotically approaches a 
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split-monopole form at far distances. Their results also show that the solution is 
smooth across the light cylinder, and there is no evidence of a singularity of the flow 
and field. In a subsequent paper™ Contopoulos and Kazanas proposed a possible 
conversion mechanism of pulsar power into particle energy in the wind region, using 
their solution. 

For the above reason, the CKF solution is quite interesting. In particular, the 
non-singular property is in stark contrast with the previous understanding that a 
dissipation-free flow would lead to a singularity in the flow and the magnetic field, 
occurring a short distance beyond the light cylinder (e.g., Ref. QJ). The drift ap- 
proximation is expected to break down near the light cylinder (e.g., Ref. ITU])) . since 
the drift velocity \v±\ = c\E x B\/B 2 approaches the light speed c due to the fact 
that \E\ ~ \B\ in the wind region. These arguments sound reasonable, but depend 
strongly on the electro-magnetic field configuration. If there were no drawbacks of 
the CKF solution, it would allow for the understanding of the pulsar magnetosphere 
to progress significantly, or it would reveal the need to change the current picture 
presented in the literaturePJ^J On the other hand, if the validity of the CKF so- 
lution could be checked, this solution could be used as the lowest approximation, 
and various physical effects could be added to it. Furthermore, it would help to 
understand other astrophysical phenomena, because the pulsar equation can also be 
applied to the black hole magnetosphere and jets (e.g., Ref. 14)). Therefore, it is 
quite important to determine the properties of the CKF solution. 

The purpose of this paper is to explore the physical quantities described by the 
CKF solution. In particular, we explicitly calculate the electro-magnetic structure 
and drift velocity. These quantities were not given in Ref. 15). We calculated 
them by differentiating the magnetic flux function. For the numerical differentiation, 
it was necessary to construct a flux function using many more grid points than 
in the computation reported in Ref. 15). Furthermore, the validity of the CKF 
solution is investigated. In §2, we briefly summarize the basic equations for the 
axially symmetric MHD field in the force-free limit. In §3, we briefly summarize the 
numerical method and the CKF solution constructed using numerical calculations. 
In §4, the magnetic field structure and drift velocity are calculated. Finally, some 
remarks are given in §5. In the appendix, we also briefly summarize some properties 
of the split-monopole magnetic field, which can be calculated analytically, and is 
useful for comparison. 



§2. Assumptions and formalism 



In this section we briefly summarize the equations governing the stationary ax- 
isymmetric magnetosphere. In the axisymmetric case, it is convenient to decouple 
the poloidal and toroidal parts of the magnetic fields in cylindrical coordinates, 
(R, <p, Z). The poloidal field B p = (Br, Bz) is written in terms of the flux function 
<P(R, Z), while the toroidal field is given by another function, A(R, Z), as 

(B R , B z ) = (~dz$, -5 A ' t 2 ' 1 ) 
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On the basis of some assumptions, the function A should depend only on ^, as dis- 
cussed below. The plasma is assumed to be frozen to the magnetic field everywhere. 
In the limit of large magnetic Reynolds number, the velocity component perpen- 
dicular to the magnetic field is described by the so-called " E x B " drift velocity, 

cE x B , , 

v± = -p-. (2-2) 

Including the velocity along the magnetic field line, the electric and magnetic fields 
satisfy the relation 

E + - x B = 0. (2-3) 

c 

This condition holds if there exists a sufficient amount of charge density, whose origin 
in discussed in the literatureP^'^ 

As a first step to understanding the overall structure without a realistic dissi- 
pation mechanism, we here ignore the small gap regions and consider the electro- 
magnetic structure on a much larger length scale. From the ideal MHD condition 
(12-31) . the electric field can be expressed by 

(E R , E+, E Z ) = (-^d R $, 0, - ^dz$j = -V (^$J , (2-4) 

where the angular velocity of the magnetic field line Q is also a function of <P in 
general, but we assume it to be constant in this paper. Therefore, the magnetosphere 
rotates like a rigid body. 

The charge density p e is given by the divergence of the electric field (Gauss's 
law), and the current j by the rotation of the magnetic field (Ampere's law). They 
are explicitly written as 

ft = -4^* = -4^(l W * # + a K)' <2 ' 5> 
0* U. h) = £ - i - |«*) , ia«A) . (2-6) 

We assume that the electro-magnetic force dominates. That is, we ignore the 
gravity, pressure and the inertial terms. The force-free condition is written 

Pe E+-xB = 0. (2-7) 
c 

After some algebra, we find that the function A is on a surface of constant <P, i.e. 
A = A(<P). The function <P should satisfy the so-called 'pulsar equation', 

(1 - f) A$ - ^d R <P + AA' = 0, (2-8) 

where [3 is the cylindrical radius normalized with respect to the light cylinder: f3 = 
RQ/c = R/Rlc- Using Eq. (|2-8|) . the current and charge density can be rewritten 
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as 



c / . 1 2B 2 B Z - AA! tl , 
Ur, U, Jz) = ^ U'B R , - - , A'B Z ) , (2-9) 



Q 2B Z -AA' 
Pe ~~4^c 1-/3 2 ' 

The current flow within the flux tube ^ is given by 



(2-10) 



where we have defined <I? = as the Z-axis. Therefore, 

A{0) = (2-12) 
is obtained from the regularity of the toroidal magnetic field (see • 

§3. Numerical scheme to solve the singular equation 

3.1. Boundary conditions 

The pulsar equation (|2-8|) is seen to have a singular surface, i.e. the light cylinder 
R = Rlc = c/f2. Both inside and outside the light cylinder, the equation becomes a 
partial differential equation of the elliptic type for a given function A. The procedure 
for constructing a numerical solution of (|2-8[) is to solve the equation on each side 
and to match the solutions at the light cylinder. The interior and exterior solutions 
are not in general continuous across the light cylinder for an assumed function A. 
Therefore, the function A should be modified iteratively so as to converge. 

We set boundary conditions for the elliptic equation that are almost the same as 
for the CKF solution. The function ^> near the star is described by a dipole magnetic 
field with magnetic dipole moment m as 

mR 2 

{R 2 + Z 2 f/ 2 ' (3 ' 1} 

The last open field line for this dipole field corresponds to <P = <P pc = m/R^c, which 
passes the light cylinder on the equatorial plane. At the polar axis, is given by 
the Dirichlet condition as 

<5(0,Z) = 0. (3-2) 

The boundary condition at the outer right points is given by the Neumann condition 

as 

d R $(R max , Z) = 0, (3-3) 
while, the upper outermost boundary is given by 

$(R, Z max ) = 0, (3-4) 
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where the domain of calculation is limited to < R < R max , < Z < Z max . 
The conditions 2? max = oo and Z max = oo, were applied to the CKF solution, 
using a certain coordinate transformation, but we limited the spatial domain as 
-Rmax = Zmax = 9Rlc- The effect of this change was numerically checked, and we 
found that there is no substantial difference between the results obtained in the two 
cases. Our choice is better to facilitate numerical convergence. We also changed 
the outer boundary conditions, but we found that they are not so sensitive to the 
boundary conditions, and the calculations converge to the result almost uniquely. 
The boundary condition on the equatorial plane is divided by the light cylinder. 
Because particles can rotate with the star within the light cylinder, we have 

d z <P(R,0)=0. (R<R LC ) (3-5) 

Outside the light cylinder, the line is extended to infinity. We therefore set it as 

#(22,0) =# p en , (R>Rlc) (3-6) 

where # pen is the last inner open field line, i.e. <&{Rlc-, 0) calculated for the interior 
solution. 

In order to avoid the singular surface of Eq. ()2-8|) , we require at the light cylinder 
R = Rlc = c/ f2 the relation 

AA' = = 2B Z . (3-7) 

In this case, the singularities of the charge density p e and toroidal current vanish 
simultaneously [See Eqs. JESHEHDJl-] 

It is not clear a priori which functional form of A(<P) satisfies all of the boundary 
conditions. We therefore need to carry out an iterative process to obtain the solution. 
Contopoulos et al. 15 ' successfully solved the equation by searching for a suitable 
function A so as to match the interior and exterior solutions at the surface. We used 
the same relaxation-type technique. 

3.2. Results 

We solved Eq. (|2-8|) with the boundary conditions in the numerical domain 
2?max = ^max = QRlc- A self-consistent solution for <P and A(<P) was calculated on a 
grid with 80 x 80 points inside and another 80 x 80 points outside the light cylinder 
to insure numerical accuracy. The calculation required more than 10 4 iterations 
to realize a smooth result. We performed many trials to search for different types 
of solutions, but just one, the CKF solution, was found. This may be a unique 
solution satisfying the boundary conditions, although we have not yet proved this. 
The overall structure of the solution is the same as that of the CKF solution, but 
the value of the last open field line is a little bit larger than that of the actual CKF 
result. The result in that case, obtained with a 30 x 30 grid, is # pen = 1.36# pc , 
while our result is # pcn = 1.66# pc . The origin of this discrepancy is presumably the 
difference in numerical resolution. (We used more than twice as many grid points 
in each direction.) We also checked their result with a calculation using a 30 x 
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30 grid. We found that the value ^ pen increases slightly with the grid points, but 
our numerical result seems to converge to a value around ^ ope n = 1.66# pc . It is 
interesting to compare this value with that for the case of constant ^4, 10 ' in which 
^open = 1.592# pc . The precise value depends on the numerical accuracy, and scales 
the physical quantities of the solution. The calculated function <P(R, Z) is shown 
in Fig. 1. Both interior and exterior functions connect smoothly across the light 
cylinder. We only show the result for the inner part of the numerical calculation, 
where < R/Rlc < 4, < Z/Rlc < 3. The result in this region is almost 
unaffected by our choice of the outermost boundary conditions at i? max , Z m£H . 

The function needed during the numerical iterations is not A, but A A' . The 
conditions adopted in our numerical calculation are A A' = for <P = and ^ > # opeil , 
where we used Eq. (|2-12j) and the implicit assumption that the electric circuit is closed 
in the whole system: 



A($)=0. (<2>><£ pen) (3-8) 



The function A is constructed from the numerical data for A A 1 ' . The result is shown 
in Fig. 2. The result indeed shows that ^4(^ open ) ^ 0. This causes serious problems 
when applying our computation to physical situations. The current closure is not 
realized, and there is a discontinuity of in this model. The total current from 
the upper hemisphere has a definite sign, J < 0. A simple way to use this model 
is to introduce a thin current sheet on ^ pen- This artificially ensures that a return 
current sheet equal to — ^4(^ open ) flows along the equatorial plane (R > Rlc) an d the 
last open field line. The discontinuity of on ^ pen is clear, because B^ / along 
the open field line, and B^ = A(<P)/R = on the closed line. This discontinuity 
may be related to the current sheet. However, in a model with the current sheet, an 
enormously strong force, such as that due to a pressure, would be required. 

One may think that a stringent condition, such as A = at <P = and ^ pcn = 0, 
should be imposed in order to construct a more realistic model (without an artificial 
current sheet). We used this condition in the numerical iteration, but we did not 
obtain a convergent solution. Our numerical experiment may not be complete, but 
it suggests that a global solution of the pulsar equation can be obtained only under 
the weak conditions on A used in previous worksP^ 

We also compared the functional form of our result with that in the split- 
monopole case. The interesting point regarding the CKF solution is the existence 
of separated regions with different signs of A'. The poloidal current flows along a 
constant <P surface, and the current in the interval d<l> is dJ = cA'd<3>/2. The region 
in which A' < (0 < < 1.36<^ pc ) corresponds to an out-flowing electron (or an 
in-flowing positron), and that in which A' > (1.36^ pc <<£> < ^ pen) corresponds to 
an in-flowing electron. However, the amount of flow in the region satisfying A' > 
in this model is insufficient to satisfy the current closure condition, J = 0. 



Properties of Axisymmetric Pulsar Magnetosphere 



7 



3 




2 



M 



1 











1 



2 
R 



3 



4 



Fig. 1. Structure of the axisymmetric force- free magnetosphere of an aligned rotating magnetic 
dipole. The solid curves represent the flux surfaces in intervals of 0.2<P pc , with $ = along 
the axis, while dashed curves represent those of the split-monopole solution. The dotted curve 
represents the last open field line, <? op cn = 1.66<P pc . The boundary condition near the star is 
imposed at the quarter circle indicated in the lower left-hand corner of the figure. The dotted- 
dash curve indicates the null surface. 



4.1. Global magnetic field structure 

In Fig. 3, we display the magnetic fields derived from the numerical solution. 
The field strength of each component is represented by the contours. We also show 
the magnetic fields for the split-monopole solution for comparison. One obviously 
different region is inside the last closed line. The overall magnetic field structure 
is almost the same, except in certain regions. This seems reasonable, since the 
difference between the fields in the functions <P and A is not so large, as shown in §3. 

We now consider the differences in detail. For Br of the CKF solution, there 
exists a discontinuity on the equatorial plane around (Rlc, 0). This behavior results 
from the boundary condition on the equatorial plane. As described in §3, the nu- 
merical boundary conditions are divided by the point {Rlc, 0)- The condition inside 
the light cylinder ensures Br = on the equatorial plane and Br / in general 
outside the light cylinder. Therefore, it might be a subtle problem to devise a proper 
numerical treatment around the point (Rlc,0). 

Constructing numerical treatment around the Z-axis may also be subtle. In the 
numerical calculation, we imposed the condition A — > as <P — > (i.e. as R —> 0), 
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Fig. 2. The electric current distribution along open field lines. The solid curve represents our 
numerical results. It is similar to the CKF solution, but is stretched to larger values of $ pc . 
The dotted curve represents the results for the split-monopole solution. 



but we have to determine the value = A/R for R s=s 0. The regular field 
should approach 0. In our numerical calculation, we do not impose such a strong 
condition as A/R — > 0. Therefore, there may be cusps on the zT-axis. Similarly, the 
Bz field may be singular on the zT-axis. We set the boundary condition as <P = 
on the axis and numerically confirmed that 8r^> ~ 0, but there is no need for the 
smoothness of the numerical value Bz = dji<P/R near R = 0. Our numerical results 
show that there is possibly singular behavior of physical quantities, Br near Z = 0, 
and Bfy and Bz near R w 0, but we need more precise calculations to reach a definite 
conclusion. 

In Fig. 4, we show the ratio of the toroidal to poloidal components of the 
magnetic field, B§/{B^ + i?^) 1 / 2 = Bt/B p . Regarding the split-monopole field, 
the toroidal component gradually increases with R and dominates completely for 
R > Rlc- It also can be seen from our result that the toroidal field dominates, but 
not as strongly, in the exterior of the light cylinder. 
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Fig. 3. Contours of the magnetic field strength for each component of the split-monopole field 
(a)-(c) and our numerical results (d)-(f). (a) The solid curves correspond to Br = 7.0, 3.5, 2.0, 
1.0, 0.5, 0.3 and 0.1B o = m/R^c- (b) The solid curves correspond to B^ = —3.5, —3.0, —2.5, 
-2.0, -1.5, -1.0 and -0.5B o . (c) The solid curves correspond to B z = 1.7, 1.5, 1.3, 1.1, 0.9, 
0.7, 0.5, 0.3 and 0.1B o . (d) The solid curves correspond to B R = 7.0, 3.5, 2.0, 1.0, 0.5, 0.3 and 
0.1B o . (e) The solid curves correspond to B^ = -2.5, -2.0, -1.5, -1.0 and -0.5B o . (f) The 
solid curves correspond to Bz = 1.3, 1.1, 0.9, 0.7, 0.5, 0.3 and 0.1B o , while dotted curves show 
outwardly Bz = —0.7, —0.5, —0.3 and —0.1 B . The null line, along which Bz = 0, is indicated 
by the dotted-dash curve. 



4.2. Drift velocity 

The components of the drift velocity ()2-2|) are explicitly given by 



*±=(-^B R B„ Mpl + Bl), -2°B z B+y (4-1) 
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Fig. 4. Contours of the ratio of the toroidal to poloidal components of the magnetic field. The 
split-monopole solution is shown in the left panel and our result is shown in the right panel. 
The contour level increases from Bt/B p — 0.5 (innermost) at intervals of 0.5. The magnitude 
of the toroidal component exceeds that of the poloidal component beyond the light cylinder in 
both solutions. 



The absolute magnitude is 



It is clear that the condition v ± < c holds inside the light cylinder. However, this 
condition is not guaranteed outside the light cylinder, and whether or not it holds de- 
pends significantly on the magnetic field configuration. Indeed, this condition holds 
everywhere for the magnetic field of the split-monopole, as shown in the appendix. 
In Fig. 5, we show the drift velocity v± calculated from the CKF solution. The nu- 
merical result shows that v± < 0.8c within the light cylinder and that the value v± 
gradually increases with R. The drift velocity exceeds the light speed for R > 3Rlc- 
This violation is clear, especially for small Z. The condition v± < c means that 
the toroidal mag netic field should dominate for B\j{B 2 R + B%) > (R/Rlc) 2 - 1, 
from Eq. I|4-2|) . The ratio of the toroidal to poloidal components, B^j(B 2 R + B 2 Z ) 1 / 2 , 
increases with R, but the dependence is not as steep as for the split-monopole field. 
This difference leads to the violation of causality, i.e. v j_ > c in the CKF solution, 
but no violation in the split-monopole case. It is thus seen that the applicable range 
of the CKF solution is limited within a few times the size of the light cylinder. 

One may expect the breakdown of the drift approximation, as suggested in the 
literaturepl but it should be noted that the breakdown depends significantly on 
the adopted electro-magnetic fields. A counterexample is the split-monopole case. 
In this case, the drift approximation holds everywhere. The breakdown of the drift 
approximation is likely to occur outside the light cylinder for most electro-magnetic 
fields. In general, the location of the breakdown is not known a priori. The drift 
velocity cannot exceed the light velocity in a realistic situation. The inertial term 
gradually prevents the increase of the drift velocity as a particle moves away from a 
star. The magnetic field near a star is strong enough that the inertial term of the 
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Fig. 5. Contours of the drift velocity for the split-monopole solution and our numerical results. 
The solid lines in the left panel indicate \v±\ = 0.8, 0.9, 0.95 and 0.97 c. The drift velocity of 
our result is shown in right-hand panel for the cases \v± \ = 0.8, 0.9, 1.1 and 1.2c (solid), and 
= 1.0c (dotted). 



particle can be ignored, but it decreases with the distance from the star. For this 
reason, the inertial term can no longer be ignored beyond a certain distance. This is 
a general picture, but the position beyond which the inertial term cannot be ignored 
depends strongly on the magnetic field configuration. In this paper, we have for the 
first time elucidated the properties of electro-magnetic fields described by the CKF 
solution and the drawbacks of the solution. 

§5. Remarks 

In this paper, we have explored the axially symmetric magnetosphere around an 
aligned rotating dipole field. The ideal MHD and force-free conditions are satisfied 
everywhere by assuming a sufficient charge density. The electro-magnetic structure 
can be determined using the so-called pulsar equation, which is a coupled equation 
for <P and A. Some authors have examined some special cases in the past, e.g. that 
in which A is constant."-' linear in (^JEa and quadratic in <&f* However, these do not 
correspond to realistic situations. A new numerical solution reported by Contopoulos 
et alP^ seems to be promising. That solution may be useful as a starting point for 
a more realistic pulsar magnetosphere. From this point of view, we have examined 
the magnetic field and drift velocity in detail. 

The overall global structure of the magnetosphere we have studied, is quite 
similar to that of a split-monopole field. Some difference in the current distribution 
A appears in regions corresponding to the case in which the field lines extend above 
the equatorial plane. A modification is necessary to connect the inner dipole-like field 
with the outer monopole-like field. This causes return currents, although the total 
amount of such currents is insufficient to realize current closure. For this reason a 
current sheet is required in the magnetosphere. Our result also shows the breakdown 
of the drift approximation beyond a few times the size of the light cylinder, and 
this result may modify the understanding regarding the CKF solution proposed in 
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Ref. I16|) . In that work, it was shown that the Lorentz factor of the outflowing 
plasma increases linearly with the distance from a neutron star. Therefore, at a 
large distance, R ~ W 3 Rlc, the inertial term can no longer be ignored (i.e., the 
force- free condition breaks down), whereas our numerical calculation shows that the 
validity of the CKF solution is limited within a short distance ~ 2>Rlc- Therefore 
some assumptions used in this calculation should incorporate the effects of a non- 
ideal MHD or the inertia of the particles in order to obtain a more realistic solution 
(e.g., see Ref.lU). 
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Appendix: Split-monopole solution 

In this appendix, we consider the physical structure represented by the split- 
monopole solution. It may be useful to compare this explanation with the CKF 
solution. As noted by CKF, their solution asymptotically approaches the split- 
monopole solution. 

When the function A is quadric in ^, the resultant magnetic field corresponds to 
the magnetic monopole solution. The solution of Eq. (|2-8j) can be explicitly written 

as 



A = -—\2- -J, (5-2) 

where q is related to the monopole charge. This fact is easily checked by calculating 
the magnetic field given as 

(R H P 1 - ( * R ggg * Z \ 

{ R > z) V0R 2 + z 2 ) 3/2 ' c(R 2 + z 2 y (R^ + z 2 f/ 2 )' 1 j 

The poloidal component of magnetic field (Br, Bz) is always radial. This poloidal 
field dominates near the origin. By contrast, the toroidal field dominates across the 
light cylinder; that is, B$ > B p = (B\ + B^) 1 / 2 for R > c/Q. In this way, this 
monopole solution represents a smooth transition of the magnetic field through the 
light cylinder. It may be helpful in obtaining physical understanding to calculate 
the current flow, 

RZ Z 2 \ 

(3R, 3<t>, Jz) - ^ ^(fl2 + Z 2)2> °> (jp + Z 2y ) ■ ( 54 ) 

The current flow is purely poloidal and out-flowing, i.e. J > 0. This current grad- 
ually enhances the toroidal magnetic field B^, which dominates outside the light 
cylinder. 
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The drift velocity is given by 

c(3 ( Rf3 ZP \ 

{V R , V(t> ,Vz) ( 1+/ g2) ^ (jR 2 + Z 2)l/2' X ' (^2 + ^2)1/2^ ■ ^ 

It is easily verified that uj_ < 1 everywhere. This drift velocity asymptotically 
approaches the light speed as R — » oo. 

On the basis of the above consideration, there seems to be no problem with the 
monopole solution, except that it may not exist in nature. No monopole has yet been 
discovered in nature, and so this solution is used only for, e.g., upper hemisphere 
with +q and the lower hemisphere with —q. Combining these two solutions with 
opposite charge, the so-called split-monopole solution is sometimes used for physical 
analysis, e.g., in Ref. 110)1 . In that case, the particle velocity of consistent flow is c 
everywhere. 
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